The paper is an addition to the growing body of research focused on the relationship between ancient societies and rapid climatic events, specifically the 4.2 kya event (a period of increased aridity and cooling). Lawrence et al. (2021) aim to take an empirical approach at a large spatial scale to investigate trends in population and settlement organization related to the impact of the 4.2 kya event across the entirety of the Northern Fertile Crescent from 6,000 to 3,000 cal. BP.
The paper employs three different datasets (archaeological data and radiocarbon dates) from the period between 6000 and 3000 cal. BP to investigate population trends. There are two different archaeological settlement datasets used in this study. The first consists of data collected from 16 published archaeological surveys and includes 1157 sites divided into 2783 occupation phases. The second consists of urban settlement data (defined as sites which reach a size of ten hectors or more) and includes 132 urban sites divided into 283 occupation phases. The third dataset is comprised of 963 radiocarbon dates from archaeological contexts collected over a slightly broader time range, 6500 to 2500 uncal. BP, in the aim to avoid edge effects.
The archaeological settlement data (both sets) were binned into a series of 100-year time slices starting at 6,000 BP and ending at 3,000 BP. The authors applied aoristic approaches to deal with uncertainties in the dataset. Essentially, this method assumes that the total probability of an archaeological event within a given time span is 1, which indicates absolute certainty. This is then divided by the length of the sites chronological range to represent the probability of existence for each temporal block (100 years) and an aoristic weight is produced (displayed in Fig. 4 and 5). In addition, Monte Carlo methods were applied to generate randomised start of occupation periods that could be compared with other demographic proxies.
The radiocarbon dataset was processed using the R package rcarbon. The dates were calibrated and a summed probability distribution (SPD) of individual dates was produced. A KDE map of the spatial distributions of radiocarbon dates was created. they compared the summed normalised probability distributions and observed unnormalized SPDs of calibrated radiocarbon dates against a theoretical null model. They fitted a logistic growth model to the observed SPD and produced a 95% confidence envelope of 1,000 SPDs to test if the observed pattern differed from the null models (displayed in Fig. 2).
The conclusions of this study indicate support for claims of a decline in population and social complexity in relation to the 4.2 kya event. There are obvious negative trends in all three datasets between 4,300 and 3,900 cal. BP and sharp drops in the SPD and rural settlement proxies.
This report will replicate almost all parts of the analyses conducted in the original paper on all three assemblages. However, it will re-frame from replicating the spatial data maps.
The archaeological datasets:
Both archaeological datasets were analysed using the same methodology. The settlement data was loaded in and the dates associated with sites were converted from BC/AD to BP for analysis and site duration was also calculated and added to the dataset. The paper split the data into two groups, above or below 300mm annual average rainfall, and so this was replicated in the analyses. Bins of 100 year time slices were created. Aoristic weights of all sites were calculated and placed into a new data frame, this was replicated with sites above 300mm or below 300mm of rainfall. The data was sumarised in order for it to be plotted, this included extracting the site counts and aoristic weights for both rainfall groups and extracting the dates of value (between 3000 - 6000 BP) and converting them into vectors that could be plotted. - Fig 4 and 5 were reproduced.
The radiocarbon dataset:
The radiocarbon data will be processed using the ‘rcarbon’ R package outline in the paper (and above). The goal of this will be to replicate the calibration and summed probability distributions (SPD) in order to reproduce Figure 2b, c, and d. The dataset will be divided into two groups. The first will be all dates and the second will be the dates produced from short-lived organic samples (e.g. bones, collagen, seeds, etc). These groups will be calibrated using the ‘calibrate()’ function and, as discussed in the paper, they will be un-normalized to avoid artificial peaks produced by normalization. SPDs of all individual dates within a time range of 6,500 to 2,500 BP are produced next. Bins of 50 years were created in order to reduced potential biases associated with SPDs. Distributions of both groups (all dates and short-lived) were produced, as well as normalized and un-normalized SPDs for comparison. Additionally, Pearsons coefficients were also calculated and correlated to 0.96 as per the paper.
sites <- read.csv("settlements.csv", sep = ",", header = T) #load in site data
head(sites) #check first lines
## Name Id StartDate EndDate Longitude Latitude Rain_mm
## 1 BON_1 2500 2100 36.85311 34.75535 318.291
## 2 BON_1 1900 1600 36.85311 34.75535 318.291
## 3 BON_1 1600 1200 36.85311 34.75535 318.291
## 4 BON_1 1000 700 36.85311 34.75535 318.291
## 5 BON_11 2500 2100 36.84555 34.77783 318.291
## 6 BON_11 1900 1600 36.84555 34.77783 318.291
## Source
## 1 Morandi Bonacossi 2007
## 2 Morandi Bonacossi 2007
## 3 Morandi Bonacossi 2007
## 4 Morandi Bonacossi 2007
## 5 Morandi Bonacossi 2007
## 6 Morandi Bonacossi 2007
nrow(sites) #2783 occupation phases (corresponds with paper)
## [1] 2783
#dates are in BC/AD and need to be converted to BP for analysis
StartdatesBP <- (sites$StartDate) + 1950 #create new variable
EnddatesBP <- (sites$EndDate) + 1950
sites[, "StartdatesBP"] <- StartdatesBP #add variables to df
sites[, "EnddatesBP"] <- EnddatesBP
#duration of sites
Duration <- sites$StartdatesBP - sites$EnddatesBP
sites[, "Duration"] <- Duration #add to df
head(sites)
## Name Id StartDate EndDate Longitude Latitude Rain_mm
## 1 BON_1 2500 2100 36.85311 34.75535 318.291
## 2 BON_1 1900 1600 36.85311 34.75535 318.291
## 3 BON_1 1600 1200 36.85311 34.75535 318.291
## 4 BON_1 1000 700 36.85311 34.75535 318.291
## 5 BON_11 2500 2100 36.84555 34.77783 318.291
## 6 BON_11 1900 1600 36.84555 34.77783 318.291
## Source StartdatesBP EnddatesBP Duration
## 1 Morandi Bonacossi 2007 4450 4050 400
## 2 Morandi Bonacossi 2007 3850 3550 300
## 3 Morandi Bonacossi 2007 3550 3150 400
## 4 Morandi Bonacossi 2007 2950 2650 300
## 5 Morandi Bonacossi 2007 4450 4050 400
## 6 Morandi Bonacossi 2007 3850 3550 300
#the paper splits the date into two groups, above or below 300mm rainfall avg.
above300mm <- subset(sites, sites$Rain_mm > 300) #variable will be used later
below300mm <- subset(sites, sites$Rain_mm <= 300)
#basic descriptive stats of the rainfall groups
summary(above300mm$Rain_mm) #rainfall stats
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 302.9 431.8 464.8 488.1 548.9 755.5
summary(below300mm$Rain_mm)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -9999.0 250.9 258.5 186.8 286.3 299.1
summary(above300mm$Duration) #duration of site occupation in both groups
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 100.0 300.0 500.0 697.1 800.0 3800.0
summary(below300mm$Duration)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 100.0 400.0 500.0 542.8 650.0 3800.0
The settlement data is binned into a series of 100-year time slices and a new column is created to put the binned results in.
binwidth <- 100
breaks <- seq(1900, 8800, by = binwidth) #include range of dates
mids <- seq(min(breaks) + (binwidth/2), max(breaks) - (binwidth/2), by=binwidth)
newcolumns <- paste(breaks[1:(length(breaks)-1)], "start", sep = "")
An aoristic method is applied to the dataset. First a new dataframe needs to be created for the results, this includes a new column for site counts and one for the aoristic weights. Then a loop is used to create aoristic weights. This process assumes that the probability of site occupation occurring at each site is 1, then divide this by the length of the sites occupation to produce a ‘weight’ which is then fed into the df created.
This process is then repeated for sites above and below 300mm of annual rainfall.
#aoristic weights for all sites
#create df for results
aoristicweights <- sitecounts <- sites
sitecounts[newcolumns] <- NA
sitecounts<-sitecounts[,newcolumns]
aoristicweights[newcolumns]<-NA
aoristicweights<-aoristicweights[,newcolumns]
aoristicweights<- sitecounts
for (i in 1:nrow(sites)){
cat(paste(i,"; ",sep=""))
sitestart <- sites$StartdatesBP[i]
siteend <- sites$EnddatesBP[i]
siteyears <- seq(siteend, sitestart, by = 1)
siteyearshist <- hist(siteyears, breaks=breaks, plot = F)
timeweights <- siteyearshist$counts/binwidth
sitecounts[i,] <- timeweights #counts
aoristicweights[i,] <- timeweights/sum(timeweights) #aoristic weights
}
## 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80; 81; 82; 83; 84; 85; 86; 87; 88; 89; 90; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100; 101; 102; 103; 104; 105; 106; 107; 108; 109; 110; 111; 112; 113; 114; 115; 116; 117; 118; 119; 120; 121; 122; 123; 124; 125; 126; 127; 128; 129; 130; 131; 132; 133; 134; 135; 136; 137; 138; 139; 140; 141; 142; 143; 144; 145; 146; 147; 148; 149; 150; 151; 152; 153; 154; 155; 156; 157; 158; 159; 160; 161; 162; 163; 164; 165; 166; 167; 168; 169; 170; 171; 172; 173; 174; 175; 176; 177; 178; 179; 180; 181; 182; 183; 184; 185; 186; 187; 188; 189; 190; 191; 192; 193; 194; 195; 196; 197; 198; 199; 200; 201; 202; 203; 204; 205; 206; 207; 208; 209; 210; 211; 212; 213; 214; 215; 216; 217; 218; 219; 220; 221; 222; 223; 224; 225; 226; 227; 228; 229; 230; 231; 232; 233; 234; 235; 236; 237; 238; 239; 240; 241; 242; 243; 244; 245; 246; 247; 248; 249; 250; 251; 252; 253; 254; 255; 256; 257; 258; 259; 260; 261; 262; 263; 264; 265; 266; 267; 268; 269; 270; 271; 272; 273; 274; 275; 276; 277; 278; 279; 280; 281; 282; 283; 284; 285; 286; 287; 288; 289; 290; 291; 292; 293; 294; 295; 296; 297; 298; 299; 300; 301; 302; 303; 304; 305; 306; 307; 308; 309; 310; 311; 312; 313; 314; 315; 316; 317; 318; 319; 320; 321; 322; 323; 324; 325; 326; 327; 328; 329; 330; 331; 332; 333; 334; 335; 336; 337; 338; 339; 340; 341; 342; 343; 344; 345; 346; 347; 348; 349; 350; 351; 352; 353; 354; 355; 356; 357; 358; 359; 360; 361; 362; 363; 364; 365; 366; 367; 368; 369; 370; 371; 372; 373; 374; 375; 376; 377; 378; 379; 380; 381; 382; 383; 384; 385; 386; 387; 388; 389; 390; 391; 392; 393; 394; 395; 396; 397; 398; 399; 400; 401; 402; 403; 404; 405; 406; 407; 408; 409; 410; 411; 412; 413; 414; 415; 416; 417; 418; 419; 420; 421; 422; 423; 424; 425; 426; 427; 428; 429; 430; 431; 432; 433; 434; 435; 436; 437; 438; 439; 440; 441; 442; 443; 444; 445; 446; 447; 448; 449; 450; 451; 452; 453; 454; 455; 456; 457; 458; 459; 460; 461; 462; 463; 464; 465; 466; 467; 468; 469; 470; 471; 472; 473; 474; 475; 476; 477; 478; 479; 480; 481; 482; 483; 484; 485; 486; 487; 488; 489; 490; 491; 492; 493; 494; 495; 496; 497; 498; 499; 500; 501; 502; 503; 504; 505; 506; 507; 508; 509; 510; 511; 512; 513; 514; 515; 516; 517; 518; 519; 520; 521; 522; 523; 524; 525; 526; 527; 528; 529; 530; 531; 532; 533; 534; 535; 536; 537; 538; 539; 540; 541; 542; 543; 544; 545; 546; 547; 548; 549; 550; 551; 552; 553; 554; 555; 556; 557; 558; 559; 560; 561; 562; 563; 564; 565; 566; 567; 568; 569; 570; 571; 572; 573; 574; 575; 576; 577; 578; 579; 580; 581; 582; 583; 584; 585; 586; 587; 588; 589; 590; 591; 592; 593; 594; 595; 596; 597; 598; 599; 600; 601; 602; 603; 604; 605; 606; 607; 608; 609; 610; 611; 612; 613; 614; 615; 616; 617; 618; 619; 620; 621; 622; 623; 624; 625; 626; 627; 628; 629; 630; 631; 632; 633; 634; 635; 636; 637; 638; 639; 640; 641; 642; 643; 644; 645; 646; 647; 648; 649; 650; 651; 652; 653; 654; 655; 656; 657; 658; 659; 660; 661; 662; 663; 664; 665; 666; 667; 668; 669; 670; 671; 672; 673; 674; 675; 676; 677; 678; 679; 680; 681; 682; 683; 684; 685; 686; 687; 688; 689; 690; 691; 692; 693; 694; 695; 696; 697; 698; 699; 700; 701; 702; 703; 704; 705; 706; 707; 708; 709; 710; 711; 712; 713; 714; 715; 716; 717; 718; 719; 720; 721; 722; 723; 724; 725; 726; 727; 728; 729; 730; 731; 732; 733; 734; 735; 736; 737; 738; 739; 740; 741; 742; 743; 744; 745; 746; 747; 748; 749; 750; 751; 752; 753; 754; 755; 756; 757; 758; 759; 760; 761; 762; 763; 764; 765; 766; 767; 768; 769; 770; 771; 772; 773; 774; 775; 776; 777; 778; 779; 780; 781; 782; 783; 784; 785; 786; 787; 788; 789; 790; 791; 792; 793; 794; 795; 796; 797; 798; 799; 800; 801; 802; 803; 804; 805; 806; 807; 808; 809; 810; 811; 812; 813; 814; 815; 816; 817; 818; 819; 820; 821; 822; 823; 824; 825; 826; 827; 828; 829; 830; 831; 832; 833; 834; 835; 836; 837; 838; 839; 840; 841; 842; 843; 844; 845; 846; 847; 848; 849; 850; 851; 852; 853; 854; 855; 856; 857; 858; 859; 860; 861; 862; 863; 864; 865; 866; 867; 868; 869; 870; 871; 872; 873; 874; 875; 876; 877; 878; 879; 880; 881; 882; 883; 884; 885; 886; 887; 888; 889; 890; 891; 892; 893; 894; 895; 896; 897; 898; 899; 900; 901; 902; 903; 904; 905; 906; 907; 908; 909; 910; 911; 912; 913; 914; 915; 916; 917; 918; 919; 920; 921; 922; 923; 924; 925; 926; 927; 928; 929; 930; 931; 932; 933; 934; 935; 936; 937; 938; 939; 940; 941; 942; 943; 944; 945; 946; 947; 948; 949; 950; 951; 952; 953; 954; 955; 956; 957; 958; 959; 960; 961; 962; 963; 964; 965; 966; 967; 968; 969; 970; 971; 972; 973; 974; 975; 976; 977; 978; 979; 980; 981; 982; 983; 984; 985; 986; 987; 988; 989; 990; 991; 992; 993; 994; 995; 996; 997; 998; 999; 1000; 1001; 1002; 1003; 1004; 1005; 1006; 1007; 1008; 1009; 1010; 1011; 1012; 1013; 1014; 1015; 1016; 1017; 1018; 1019; 1020; 1021; 1022; 1023; 1024; 1025; 1026; 1027; 1028; 1029; 1030; 1031; 1032; 1033; 1034; 1035; 1036; 1037; 1038; 1039; 1040; 1041; 1042; 1043; 1044; 1045; 1046; 1047; 1048; 1049; 1050; 1051; 1052; 1053; 1054; 1055; 1056; 1057; 1058; 1059; 1060; 1061; 1062; 1063; 1064; 1065; 1066; 1067; 1068; 1069; 1070; 1071; 1072; 1073; 1074; 1075; 1076; 1077; 1078; 1079; 1080; 1081; 1082; 1083; 1084; 1085; 1086; 1087; 1088; 1089; 1090; 1091; 1092; 1093; 1094; 1095; 1096; 1097; 1098; 1099; 1100; 1101; 1102; 1103; 1104; 1105; 1106; 1107; 1108; 1109; 1110; 1111; 1112; 1113; 1114; 1115; 1116; 1117; 1118; 1119; 1120; 1121; 1122; 1123; 1124; 1125; 1126; 1127; 1128; 1129; 1130; 1131; 1132; 1133; 1134; 1135; 1136; 1137; 1138; 1139; 1140; 1141; 1142; 1143; 1144; 1145; 1146; 1147; 1148; 1149; 1150; 1151; 1152; 1153; 1154; 1155; 1156; 1157; 1158; 1159; 1160; 1161; 1162; 1163; 1164; 1165; 1166; 1167; 1168; 1169; 1170; 1171; 1172; 1173; 1174; 1175; 1176; 1177; 1178; 1179; 1180; 1181; 1182; 1183; 1184; 1185; 1186; 1187; 1188; 1189; 1190; 1191; 1192; 1193; 1194; 1195; 1196; 1197; 1198; 1199; 1200; 1201; 1202; 1203; 1204; 1205; 1206; 1207; 1208; 1209; 1210; 1211; 1212; 1213; 1214; 1215; 1216; 1217; 1218; 1219; 1220; 1221; 1222; 1223; 1224; 1225; 1226; 1227; 1228; 1229; 1230; 1231; 1232; 1233; 1234; 1235; 1236; 1237; 1238; 1239; 1240; 1241; 1242; 1243; 1244; 1245; 1246; 1247; 1248; 1249; 1250; 1251; 1252; 1253; 1254; 1255; 1256; 1257; 1258; 1259; 1260; 1261; 1262; 1263; 1264; 1265; 1266; 1267; 1268; 1269; 1270; 1271; 1272; 1273; 1274; 1275; 1276; 1277; 1278; 1279; 1280; 1281; 1282; 1283; 1284; 1285; 1286; 1287; 1288; 1289; 1290; 1291; 1292; 1293; 1294; 1295; 1296; 1297; 1298; 1299; 1300; 1301; 1302; 1303; 1304; 1305; 1306; 1307; 1308; 1309; 1310; 1311; 1312; 1313; 1314; 1315; 1316; 1317; 1318; 1319; 1320; 1321; 1322; 1323; 1324; 1325; 1326; 1327; 1328; 1329; 1330; 1331; 1332; 1333; 1334; 1335; 1336; 1337; 1338; 1339; 1340; 1341; 1342; 1343; 1344; 1345; 1346; 1347; 1348; 1349; 1350; 1351; 1352; 1353; 1354; 1355; 1356; 1357; 1358; 1359; 1360; 1361; 1362; 1363; 1364; 1365; 1366; 1367; 1368; 1369; 1370; 1371; 1372; 1373; 1374; 1375; 1376; 1377; 1378; 1379; 1380; 1381; 1382; 1383; 1384; 1385; 1386; 1387; 1388; 1389; 1390; 1391; 1392; 1393; 1394; 1395; 1396; 1397; 1398; 1399; 1400; 1401; 1402; 1403; 1404; 1405; 1406; 1407; 1408; 1409; 1410; 1411; 1412; 1413; 1414; 1415; 1416; 1417; 1418; 1419; 1420; 1421; 1422; 1423; 1424; 1425; 1426; 1427; 1428; 1429; 1430; 1431; 1432; 1433; 1434; 1435; 1436; 1437; 1438; 1439; 1440; 1441; 1442; 1443; 1444; 1445; 1446; 1447; 1448; 1449; 1450; 1451; 1452; 1453; 1454; 1455; 1456; 1457; 1458; 1459; 1460; 1461; 1462; 1463; 1464; 1465; 1466; 1467; 1468; 1469; 1470; 1471; 1472; 1473; 1474; 1475; 1476; 1477; 1478; 1479; 1480; 1481; 1482; 1483; 1484; 1485; 1486; 1487; 1488; 1489; 1490; 1491; 1492; 1493; 1494; 1495; 1496; 1497; 1498; 1499; 1500; 1501; 1502; 1503; 1504; 1505; 1506; 1507; 1508; 1509; 1510; 1511; 1512; 1513; 1514; 1515; 1516; 1517; 1518; 1519; 1520; 1521; 1522; 1523; 1524; 1525; 1526; 1527; 1528; 1529; 1530; 1531; 1532; 1533; 1534; 1535; 1536; 1537; 1538; 1539; 1540; 1541; 1542; 1543; 1544; 1545; 1546; 1547; 1548; 1549; 1550; 1551; 1552; 1553; 1554; 1555; 1556; 1557; 1558; 1559; 1560; 1561; 1562; 1563; 1564; 1565; 1566; 1567; 1568; 1569; 1570; 1571; 1572; 1573; 1574; 1575; 1576; 1577; 1578; 1579; 1580; 1581; 1582; 1583; 1584; 1585; 1586; 1587; 1588; 1589; 1590; 1591; 1592; 1593; 1594; 1595; 1596; 1597; 1598; 1599; 1600; 1601; 1602; 1603; 1604; 1605; 1606; 1607; 1608; 1609; 1610; 1611; 1612; 1613; 1614; 1615; 1616; 1617; 1618; 1619; 1620; 1621; 1622; 1623; 1624; 1625; 1626; 1627; 1628; 1629; 1630; 1631; 1632; 1633; 1634; 1635; 1636; 1637; 1638; 1639; 1640; 1641; 1642; 1643; 1644; 1645; 1646; 1647; 1648; 1649; 1650; 1651; 1652; 1653; 1654; 1655; 1656; 1657; 1658; 1659; 1660; 1661; 1662; 1663; 1664; 1665; 1666; 1667; 1668; 1669; 1670; 1671; 1672; 1673; 1674; 1675; 1676; 1677; 1678; 1679; 1680; 1681; 1682; 1683; 1684; 1685; 1686; 1687; 1688; 1689; 1690; 1691; 1692; 1693; 1694; 1695; 1696; 1697; 1698; 1699; 1700; 1701; 1702; 1703; 1704; 1705; 1706; 1707; 1708; 1709; 1710; 1711; 1712; 1713; 1714; 1715; 1716; 1717; 1718; 1719; 1720; 1721; 1722; 1723; 1724; 1725; 1726; 1727; 1728; 1729; 1730; 1731; 1732; 1733; 1734; 1735; 1736; 1737; 1738; 1739; 1740; 1741; 1742; 1743; 1744; 1745; 1746; 1747; 1748; 1749; 1750; 1751; 1752; 1753; 1754; 1755; 1756; 1757; 1758; 1759; 1760; 1761; 1762; 1763; 1764; 1765; 1766; 1767; 1768; 1769; 1770; 1771; 1772; 1773; 1774; 1775; 1776; 1777; 1778; 1779; 1780; 1781; 1782; 1783; 1784; 1785; 1786; 1787; 1788; 1789; 1790; 1791; 1792; 1793; 1794; 1795; 1796; 1797; 1798; 1799; 1800; 1801; 1802; 1803; 1804; 1805; 1806; 1807; 1808; 1809; 1810; 1811; 1812; 1813; 1814; 1815; 1816; 1817; 1818; 1819; 1820; 1821; 1822; 1823; 1824; 1825; 1826; 1827; 1828; 1829; 1830; 1831; 1832; 1833; 1834; 1835; 1836; 1837; 1838; 1839; 1840; 1841; 1842; 1843; 1844; 1845; 1846; 1847; 1848; 1849; 1850; 1851; 1852; 1853; 1854; 1855; 1856; 1857; 1858; 1859; 1860; 1861; 1862; 1863; 1864; 1865; 1866; 1867; 1868; 1869; 1870; 1871; 1872; 1873; 1874; 1875; 1876; 1877; 1878; 1879; 1880; 1881; 1882; 1883; 1884; 1885; 1886; 1887; 1888; 1889; 1890; 1891; 1892; 1893; 1894; 1895; 1896; 1897; 1898; 1899; 1900; 1901; 1902; 1903; 1904; 1905; 1906; 1907; 1908; 1909; 1910; 1911; 1912; 1913; 1914; 1915; 1916; 1917; 1918; 1919; 1920; 1921; 1922; 1923; 1924; 1925; 1926; 1927; 1928; 1929; 1930; 1931; 1932; 1933; 1934; 1935; 1936; 1937; 1938; 1939; 1940; 1941; 1942; 1943; 1944; 1945; 1946; 1947; 1948; 1949; 1950; 1951; 1952; 1953; 1954; 1955; 1956; 1957; 1958; 1959; 1960; 1961; 1962; 1963; 1964; 1965; 1966; 1967; 1968; 1969; 1970; 1971; 1972; 1973; 1974; 1975; 1976; 1977; 1978; 1979; 1980; 1981; 1982; 1983; 1984; 1985; 1986; 1987; 1988; 1989; 1990; 1991; 1992; 1993; 1994; 1995; 1996; 1997; 1998; 1999; 2000; 2001; 2002; 2003; 2004; 2005; 2006; 2007; 2008; 2009; 2010; 2011; 2012; 2013; 2014; 2015; 2016; 2017; 2018; 2019; 2020; 2021; 2022; 2023; 2024; 2025; 2026; 2027; 2028; 2029; 2030; 2031; 2032; 2033; 2034; 2035; 2036; 2037; 2038; 2039; 2040; 2041; 2042; 2043; 2044; 2045; 2046; 2047; 2048; 2049; 2050; 2051; 2052; 2053; 2054; 2055; 2056; 2057; 2058; 2059; 2060; 2061; 2062; 2063; 2064; 2065; 2066; 2067; 2068; 2069; 2070; 2071; 2072; 2073; 2074; 2075; 2076; 2077; 2078; 2079; 2080; 2081; 2082; 2083; 2084; 2085; 2086; 2087; 2088; 2089; 2090; 2091; 2092; 2093; 2094; 2095; 2096; 2097; 2098; 2099; 2100; 2101; 2102; 2103; 2104; 2105; 2106; 2107; 2108; 2109; 2110; 2111; 2112; 2113; 2114; 2115; 2116; 2117; 2118; 2119; 2120; 2121; 2122; 2123; 2124; 2125; 2126; 2127; 2128; 2129; 2130; 2131; 2132; 2133; 2134; 2135; 2136; 2137; 2138; 2139; 2140; 2141; 2142; 2143; 2144; 2145; 2146; 2147; 2148; 2149; 2150; 2151; 2152; 2153; 2154; 2155; 2156; 2157; 2158; 2159; 2160; 2161; 2162; 2163; 2164; 2165; 2166; 2167; 2168; 2169; 2170; 2171; 2172; 2173; 2174; 2175; 2176; 2177; 2178; 2179; 2180; 2181; 2182; 2183; 2184; 2185; 2186; 2187; 2188; 2189; 2190; 2191; 2192; 2193; 2194; 2195; 2196; 2197; 2198; 2199; 2200; 2201; 2202; 2203; 2204; 2205; 2206; 2207; 2208; 2209; 2210; 2211; 2212; 2213; 2214; 2215; 2216; 2217; 2218; 2219; 2220; 2221; 2222; 2223; 2224; 2225; 2226; 2227; 2228; 2229; 2230; 2231; 2232; 2233; 2234; 2235; 2236; 2237; 2238; 2239; 2240; 2241; 2242; 2243; 2244; 2245; 2246; 2247; 2248; 2249; 2250; 2251; 2252; 2253; 2254; 2255; 2256; 2257; 2258; 2259; 2260; 2261; 2262; 2263; 2264; 2265; 2266; 2267; 2268; 2269; 2270; 2271; 2272; 2273; 2274; 2275; 2276; 2277; 2278; 2279; 2280; 2281; 2282; 2283; 2284; 2285; 2286; 2287; 2288; 2289; 2290; 2291; 2292; 2293; 2294; 2295; 2296; 2297; 2298; 2299; 2300; 2301; 2302; 2303; 2304; 2305; 2306; 2307; 2308; 2309; 2310; 2311; 2312; 2313; 2314; 2315; 2316; 2317; 2318; 2319; 2320; 2321; 2322; 2323; 2324; 2325; 2326; 2327; 2328; 2329; 2330; 2331; 2332; 2333; 2334; 2335; 2336; 2337; 2338; 2339; 2340; 2341; 2342; 2343; 2344; 2345; 2346; 2347; 2348; 2349; 2350; 2351; 2352; 2353; 2354; 2355; 2356; 2357; 2358; 2359; 2360; 2361; 2362; 2363; 2364; 2365; 2366; 2367; 2368; 2369; 2370; 2371; 2372; 2373; 2374; 2375; 2376; 2377; 2378; 2379; 2380; 2381; 2382; 2383; 2384; 2385; 2386; 2387; 2388; 2389; 2390; 2391; 2392; 2393; 2394; 2395; 2396; 2397; 2398; 2399; 2400; 2401; 2402; 2403; 2404; 2405; 2406; 2407; 2408; 2409; 2410; 2411; 2412; 2413; 2414; 2415; 2416; 2417; 2418; 2419; 2420; 2421; 2422; 2423; 2424; 2425; 2426; 2427; 2428; 2429; 2430; 2431; 2432; 2433; 2434; 2435; 2436; 2437; 2438; 2439; 2440; 2441; 2442; 2443; 2444; 2445; 2446; 2447; 2448; 2449; 2450; 2451; 2452; 2453; 2454; 2455; 2456; 2457; 2458; 2459; 2460; 2461; 2462; 2463; 2464; 2465; 2466; 2467; 2468; 2469; 2470; 2471; 2472; 2473; 2474; 2475; 2476; 2477; 2478; 2479; 2480; 2481; 2482; 2483; 2484; 2485; 2486; 2487; 2488; 2489; 2490; 2491; 2492; 2493; 2494; 2495; 2496; 2497; 2498; 2499; 2500; 2501; 2502; 2503; 2504; 2505; 2506; 2507; 2508; 2509; 2510; 2511; 2512; 2513; 2514; 2515; 2516; 2517; 2518; 2519; 2520; 2521; 2522; 2523; 2524; 2525; 2526; 2527; 2528; 2529; 2530; 2531; 2532; 2533; 2534; 2535; 2536; 2537; 2538; 2539; 2540; 2541; 2542; 2543; 2544; 2545; 2546; 2547; 2548; 2549; 2550; 2551; 2552; 2553; 2554; 2555; 2556; 2557; 2558; 2559; 2560; 2561; 2562; 2563; 2564; 2565; 2566; 2567; 2568; 2569; 2570; 2571; 2572; 2573; 2574; 2575; 2576; 2577; 2578; 2579; 2580; 2581; 2582; 2583; 2584; 2585; 2586; 2587; 2588; 2589; 2590; 2591; 2592; 2593; 2594; 2595; 2596; 2597; 2598; 2599; 2600; 2601; 2602; 2603; 2604; 2605; 2606; 2607; 2608; 2609; 2610; 2611; 2612; 2613; 2614; 2615; 2616; 2617; 2618; 2619; 2620; 2621; 2622; 2623; 2624; 2625; 2626; 2627; 2628; 2629; 2630; 2631; 2632; 2633; 2634; 2635; 2636; 2637; 2638; 2639; 2640; 2641; 2642; 2643; 2644; 2645; 2646; 2647; 2648; 2649; 2650; 2651; 2652; 2653; 2654; 2655; 2656; 2657; 2658; 2659; 2660; 2661; 2662; 2663; 2664; 2665; 2666; 2667; 2668; 2669; 2670; 2671; 2672; 2673; 2674; 2675; 2676; 2677; 2678; 2679; 2680; 2681; 2682; 2683; 2684; 2685; 2686; 2687; 2688; 2689; 2690; 2691; 2692; 2693; 2694; 2695; 2696; 2697; 2698; 2699; 2700; 2701; 2702; 2703; 2704; 2705; 2706; 2707; 2708; 2709; 2710; 2711; 2712; 2713; 2714; 2715; 2716; 2717; 2718; 2719; 2720; 2721; 2722; 2723; 2724; 2725; 2726; 2727; 2728; 2729; 2730; 2731; 2732; 2733; 2734; 2735; 2736; 2737; 2738; 2739; 2740; 2741; 2742; 2743; 2744; 2745; 2746; 2747; 2748; 2749; 2750; 2751; 2752; 2753; 2754; 2755; 2756; 2757; 2758; 2759; 2760; 2761; 2762; 2763; 2764; 2765; 2766; 2767; 2768; 2769; 2770; 2771; 2772; 2773; 2774; 2775; 2776; 2777; 2778; 2779; 2780; 2781; 2782; 2783;
#aoristic weights for all sites with rainfall above 300mm
aoristicweights_above300mm <- sitecounts_above300mm <- above300mm
sitecounts_above300mm[newcolumns] <- NA
sitecounts_above300mm<-sitecounts_above300mm[,newcolumns]
aoristicweights_above300mm[newcolumns]<-NA
aoristicweights_above300mm<-aoristicweights_above300mm[,newcolumns]
aoristicweights_above300mm<- sitecounts_above300mm
for (i in 1:nrow(above300mm)){
cat(paste(i,"; ",sep=""))
sitestart <- above300mm$StartdatesBP[i]
siteend <- above300mm$EnddatesBP[i]
siteyears <- seq(siteend,sitestart,by=1)
siteyearshist <- hist(siteyears, breaks=breaks, plot=FALSE)
timeweights <- siteyearshist$counts/binwidth
sitecounts_above300mm[i,] <- timeweights #counts
aoristicweights_above300mm[i,] <- timeweights/sum(timeweights) #aoristic weights
}
## 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80; 81; 82; 83; 84; 85; 86; 87; 88; 89; 90; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100; 101; 102; 103; 104; 105; 106; 107; 108; 109; 110; 111; 112; 113; 114; 115; 116; 117; 118; 119; 120; 121; 122; 123; 124; 125; 126; 127; 128; 129; 130; 131; 132; 133; 134; 135; 136; 137; 138; 139; 140; 141; 142; 143; 144; 145; 146; 147; 148; 149; 150; 151; 152; 153; 154; 155; 156; 157; 158; 159; 160; 161; 162; 163; 164; 165; 166; 167; 168; 169; 170; 171; 172; 173; 174; 175; 176; 177; 178; 179; 180; 181; 182; 183; 184; 185; 186; 187; 188; 189; 190; 191; 192; 193; 194; 195; 196; 197; 198; 199; 200; 201; 202; 203; 204; 205; 206; 207; 208; 209; 210; 211; 212; 213; 214; 215; 216; 217; 218; 219; 220; 221; 222; 223; 224; 225; 226; 227; 228; 229; 230; 231; 232; 233; 234; 235; 236; 237; 238; 239; 240; 241; 242; 243; 244; 245; 246; 247; 248; 249; 250; 251; 252; 253; 254; 255; 256; 257; 258; 259; 260; 261; 262; 263; 264; 265; 266; 267; 268; 269; 270; 271; 272; 273; 274; 275; 276; 277; 278; 279; 280; 281; 282; 283; 284; 285; 286; 287; 288; 289; 290; 291; 292; 293; 294; 295; 296; 297; 298; 299; 300; 301; 302; 303; 304; 305; 306; 307; 308; 309; 310; 311; 312; 313; 314; 315; 316; 317; 318; 319; 320; 321; 322; 323; 324; 325; 326; 327; 328; 329; 330; 331; 332; 333; 334; 335; 336; 337; 338; 339; 340; 341; 342; 343; 344; 345; 346; 347; 348; 349; 350; 351; 352; 353; 354; 355; 356; 357; 358; 359; 360; 361; 362; 363; 364; 365; 366; 367; 368; 369; 370; 371; 372; 373; 374; 375; 376; 377; 378; 379; 380; 381; 382; 383; 384; 385; 386; 387; 388; 389; 390; 391; 392; 393; 394; 395; 396; 397; 398; 399; 400; 401; 402; 403; 404; 405; 406; 407; 408; 409; 410; 411; 412; 413; 414; 415; 416; 417; 418; 419; 420; 421; 422; 423; 424; 425; 426; 427; 428; 429; 430; 431; 432; 433; 434; 435; 436; 437; 438; 439; 440; 441; 442; 443; 444; 445; 446; 447; 448; 449; 450; 451; 452; 453; 454; 455; 456; 457; 458; 459; 460; 461; 462; 463; 464; 465; 466; 467; 468; 469; 470; 471; 472; 473; 474; 475; 476; 477; 478; 479; 480; 481; 482; 483; 484; 485; 486; 487; 488; 489; 490; 491; 492; 493; 494; 495; 496; 497; 498; 499; 500; 501; 502; 503; 504; 505; 506; 507; 508; 509; 510; 511; 512; 513; 514; 515; 516; 517; 518; 519; 520; 521; 522; 523; 524; 525; 526; 527; 528; 529; 530; 531; 532; 533; 534; 535; 536; 537; 538; 539; 540; 541; 542; 543; 544; 545; 546; 547; 548; 549; 550; 551; 552; 553; 554; 555; 556; 557; 558; 559; 560; 561; 562; 563; 564; 565; 566; 567; 568; 569; 570; 571; 572; 573; 574; 575; 576; 577; 578; 579; 580; 581; 582; 583; 584; 585; 586; 587; 588; 589; 590; 591; 592; 593; 594; 595; 596; 597; 598; 599; 600; 601; 602; 603; 604; 605; 606; 607; 608; 609; 610; 611; 612; 613; 614; 615; 616; 617; 618; 619; 620; 621; 622; 623; 624; 625; 626; 627; 628; 629; 630; 631; 632; 633; 634; 635; 636; 637; 638; 639; 640; 641; 642; 643; 644; 645; 646; 647; 648; 649; 650; 651; 652; 653; 654; 655; 656; 657; 658; 659; 660; 661; 662; 663; 664; 665; 666; 667; 668; 669; 670; 671; 672; 673; 674; 675; 676; 677; 678; 679; 680; 681; 682; 683; 684; 685; 686; 687; 688; 689; 690; 691; 692; 693; 694; 695; 696; 697; 698; 699; 700; 701; 702; 703; 704; 705; 706; 707; 708; 709; 710; 711; 712; 713; 714; 715; 716; 717; 718; 719; 720; 721; 722; 723; 724; 725; 726; 727; 728; 729; 730; 731; 732; 733; 734; 735; 736; 737; 738; 739; 740; 741; 742; 743; 744; 745; 746; 747; 748; 749; 750; 751; 752; 753; 754; 755; 756; 757; 758; 759; 760; 761; 762; 763; 764; 765; 766; 767; 768; 769; 770; 771; 772; 773; 774; 775; 776; 777; 778; 779; 780; 781; 782; 783; 784; 785; 786; 787; 788; 789; 790; 791; 792; 793; 794; 795; 796; 797; 798; 799; 800; 801; 802; 803; 804; 805; 806; 807; 808; 809; 810; 811; 812; 813; 814; 815; 816; 817; 818; 819; 820; 821; 822; 823; 824; 825; 826; 827; 828; 829; 830; 831; 832; 833; 834; 835; 836; 837; 838; 839; 840; 841; 842; 843; 844; 845; 846; 847; 848; 849; 850; 851; 852; 853; 854; 855; 856; 857; 858; 859; 860; 861; 862; 863; 864; 865; 866; 867; 868; 869; 870; 871; 872; 873; 874; 875; 876; 877; 878; 879; 880; 881; 882; 883; 884; 885; 886; 887; 888; 889; 890; 891; 892; 893; 894; 895; 896; 897; 898; 899; 900; 901; 902; 903; 904; 905; 906; 907; 908; 909; 910; 911; 912; 913; 914; 915; 916; 917; 918; 919; 920; 921; 922; 923; 924; 925; 926; 927; 928; 929; 930; 931; 932; 933; 934; 935; 936; 937; 938; 939; 940; 941; 942; 943; 944; 945; 946; 947; 948; 949; 950; 951; 952; 953; 954; 955; 956; 957; 958; 959; 960; 961; 962; 963; 964; 965; 966; 967; 968; 969; 970; 971; 972; 973; 974; 975; 976; 977; 978; 979; 980; 981; 982; 983; 984; 985; 986; 987; 988; 989; 990; 991; 992; 993; 994; 995; 996; 997; 998; 999; 1000; 1001; 1002; 1003; 1004; 1005; 1006; 1007; 1008; 1009; 1010; 1011; 1012; 1013; 1014; 1015; 1016; 1017; 1018; 1019; 1020; 1021; 1022; 1023; 1024; 1025; 1026; 1027; 1028; 1029; 1030; 1031; 1032; 1033; 1034; 1035; 1036; 1037; 1038; 1039; 1040; 1041; 1042; 1043; 1044; 1045; 1046; 1047; 1048; 1049; 1050; 1051; 1052; 1053; 1054; 1055; 1056; 1057; 1058; 1059; 1060; 1061; 1062; 1063; 1064; 1065; 1066; 1067; 1068; 1069; 1070; 1071; 1072; 1073; 1074; 1075; 1076; 1077; 1078; 1079; 1080; 1081; 1082; 1083; 1084; 1085; 1086; 1087; 1088; 1089; 1090; 1091; 1092; 1093; 1094; 1095; 1096; 1097; 1098; 1099; 1100; 1101; 1102; 1103; 1104; 1105; 1106; 1107; 1108; 1109; 1110; 1111; 1112; 1113; 1114; 1115; 1116; 1117; 1118; 1119; 1120; 1121; 1122; 1123; 1124; 1125; 1126; 1127; 1128; 1129; 1130; 1131; 1132; 1133; 1134; 1135; 1136; 1137; 1138; 1139; 1140; 1141; 1142; 1143; 1144; 1145; 1146; 1147; 1148; 1149; 1150; 1151; 1152; 1153; 1154; 1155; 1156; 1157; 1158; 1159; 1160; 1161; 1162; 1163; 1164; 1165; 1166; 1167; 1168; 1169; 1170; 1171; 1172; 1173; 1174; 1175; 1176; 1177; 1178; 1179; 1180; 1181; 1182; 1183; 1184; 1185; 1186; 1187; 1188; 1189; 1190; 1191; 1192; 1193; 1194; 1195; 1196; 1197; 1198; 1199; 1200; 1201; 1202; 1203; 1204; 1205; 1206; 1207; 1208; 1209; 1210; 1211; 1212; 1213; 1214; 1215; 1216; 1217; 1218; 1219; 1220; 1221; 1222; 1223; 1224; 1225; 1226; 1227; 1228; 1229; 1230; 1231; 1232; 1233; 1234; 1235; 1236; 1237; 1238; 1239; 1240; 1241; 1242; 1243; 1244; 1245; 1246; 1247; 1248; 1249; 1250; 1251; 1252; 1253; 1254; 1255; 1256; 1257; 1258; 1259; 1260; 1261; 1262; 1263; 1264; 1265; 1266; 1267; 1268; 1269; 1270; 1271; 1272; 1273; 1274; 1275; 1276; 1277; 1278; 1279; 1280; 1281; 1282; 1283; 1284; 1285; 1286; 1287; 1288; 1289; 1290; 1291; 1292; 1293; 1294; 1295; 1296; 1297; 1298; 1299; 1300; 1301; 1302; 1303; 1304; 1305; 1306; 1307; 1308; 1309; 1310; 1311; 1312; 1313; 1314; 1315; 1316; 1317; 1318; 1319; 1320; 1321; 1322; 1323; 1324; 1325; 1326; 1327; 1328; 1329; 1330; 1331; 1332; 1333; 1334; 1335; 1336; 1337; 1338; 1339; 1340; 1341; 1342; 1343; 1344; 1345; 1346; 1347; 1348; 1349; 1350; 1351; 1352; 1353; 1354; 1355; 1356; 1357; 1358; 1359; 1360; 1361; 1362; 1363; 1364; 1365; 1366; 1367; 1368; 1369; 1370; 1371; 1372; 1373; 1374; 1375; 1376; 1377; 1378; 1379; 1380; 1381; 1382; 1383; 1384; 1385; 1386; 1387; 1388; 1389; 1390; 1391; 1392; 1393; 1394; 1395; 1396; 1397; 1398; 1399; 1400; 1401; 1402; 1403; 1404; 1405; 1406; 1407; 1408; 1409; 1410; 1411; 1412; 1413; 1414; 1415; 1416; 1417; 1418; 1419; 1420; 1421; 1422; 1423; 1424; 1425; 1426; 1427; 1428; 1429; 1430; 1431; 1432; 1433; 1434; 1435; 1436; 1437; 1438; 1439; 1440; 1441; 1442; 1443; 1444; 1445; 1446; 1447; 1448; 1449; 1450; 1451; 1452; 1453; 1454; 1455; 1456; 1457; 1458; 1459; 1460; 1461; 1462; 1463; 1464; 1465; 1466; 1467; 1468; 1469; 1470; 1471; 1472; 1473; 1474; 1475; 1476; 1477; 1478; 1479; 1480; 1481; 1482; 1483; 1484; 1485; 1486; 1487; 1488; 1489; 1490; 1491; 1492; 1493; 1494; 1495; 1496; 1497; 1498; 1499; 1500; 1501; 1502; 1503; 1504; 1505; 1506; 1507; 1508; 1509; 1510; 1511; 1512; 1513; 1514; 1515; 1516; 1517; 1518; 1519; 1520; 1521; 1522; 1523; 1524; 1525; 1526; 1527; 1528; 1529; 1530; 1531; 1532; 1533; 1534; 1535; 1536; 1537; 1538; 1539; 1540; 1541; 1542; 1543; 1544; 1545; 1546; 1547; 1548; 1549; 1550; 1551; 1552; 1553; 1554; 1555; 1556; 1557; 1558; 1559; 1560; 1561; 1562; 1563; 1564; 1565; 1566; 1567; 1568; 1569; 1570; 1571; 1572; 1573; 1574; 1575; 1576; 1577; 1578; 1579; 1580; 1581; 1582; 1583; 1584; 1585; 1586; 1587; 1588; 1589; 1590; 1591; 1592; 1593; 1594; 1595; 1596; 1597; 1598; 1599; 1600; 1601; 1602; 1603; 1604; 1605; 1606; 1607; 1608; 1609; 1610; 1611; 1612; 1613; 1614; 1615; 1616; 1617; 1618; 1619; 1620; 1621; 1622; 1623; 1624; 1625; 1626; 1627; 1628; 1629; 1630; 1631; 1632; 1633; 1634; 1635; 1636; 1637; 1638; 1639; 1640; 1641; 1642; 1643; 1644; 1645; 1646; 1647; 1648; 1649; 1650; 1651; 1652; 1653; 1654; 1655; 1656; 1657; 1658; 1659; 1660; 1661; 1662; 1663; 1664; 1665; 1666; 1667; 1668; 1669; 1670; 1671; 1672; 1673; 1674; 1675; 1676; 1677; 1678; 1679; 1680; 1681; 1682; 1683; 1684; 1685; 1686; 1687; 1688; 1689; 1690; 1691; 1692; 1693; 1694; 1695; 1696; 1697; 1698; 1699; 1700; 1701; 1702; 1703; 1704; 1705; 1706; 1707; 1708; 1709; 1710; 1711; 1712; 1713; 1714; 1715; 1716; 1717; 1718; 1719; 1720; 1721; 1722; 1723; 1724; 1725; 1726; 1727; 1728; 1729; 1730; 1731; 1732; 1733; 1734; 1735; 1736; 1737; 1738; 1739; 1740; 1741; 1742; 1743; 1744; 1745; 1746; 1747; 1748; 1749; 1750; 1751; 1752; 1753; 1754; 1755; 1756; 1757; 1758; 1759; 1760; 1761; 1762; 1763; 1764; 1765; 1766; 1767; 1768; 1769; 1770; 1771; 1772; 1773; 1774; 1775; 1776; 1777; 1778; 1779; 1780; 1781; 1782; 1783; 1784; 1785; 1786; 1787; 1788; 1789; 1790; 1791; 1792; 1793; 1794; 1795; 1796; 1797; 1798; 1799; 1800; 1801; 1802; 1803; 1804; 1805; 1806; 1807; 1808; 1809; 1810; 1811; 1812; 1813; 1814; 1815; 1816; 1817; 1818; 1819; 1820; 1821; 1822; 1823; 1824; 1825; 1826; 1827; 1828; 1829; 1830; 1831; 1832; 1833; 1834; 1835; 1836; 1837; 1838; 1839; 1840; 1841; 1842; 1843; 1844; 1845; 1846; 1847; 1848; 1849; 1850; 1851; 1852; 1853; 1854; 1855; 1856; 1857; 1858; 1859; 1860; 1861; 1862; 1863; 1864; 1865; 1866; 1867; 1868; 1869; 1870; 1871; 1872; 1873; 1874; 1875; 1876; 1877; 1878; 1879; 1880; 1881; 1882; 1883; 1884; 1885; 1886; 1887; 1888; 1889; 1890; 1891; 1892; 1893; 1894; 1895; 1896; 1897; 1898; 1899; 1900; 1901; 1902; 1903; 1904; 1905; 1906; 1907; 1908; 1909; 1910; 1911; 1912; 1913; 1914; 1915; 1916; 1917; 1918; 1919; 1920; 1921; 1922; 1923; 1924; 1925; 1926; 1927; 1928; 1929; 1930; 1931; 1932; 1933; 1934; 1935; 1936; 1937; 1938; 1939; 1940; 1941; 1942; 1943; 1944; 1945; 1946; 1947; 1948; 1949; 1950; 1951; 1952; 1953; 1954; 1955; 1956; 1957; 1958; 1959; 1960; 1961; 1962; 1963; 1964; 1965; 1966; 1967; 1968; 1969; 1970; 1971; 1972; 1973; 1974; 1975; 1976; 1977; 1978; 1979; 1980; 1981; 1982; 1983; 1984; 1985; 1986; 1987; 1988; 1989; 1990; 1991; 1992; 1993; 1994; 1995; 1996; 1997; 1998; 1999; 2000; 2001; 2002; 2003; 2004; 2005; 2006; 2007; 2008; 2009; 2010; 2011; 2012; 2013; 2014; 2015; 2016; 2017; 2018; 2019; 2020; 2021; 2022; 2023; 2024; 2025; 2026; 2027; 2028; 2029; 2030; 2031; 2032; 2033; 2034; 2035; 2036; 2037; 2038; 2039; 2040; 2041; 2042; 2043; 2044; 2045; 2046; 2047; 2048; 2049; 2050; 2051; 2052; 2053; 2054; 2055; 2056; 2057; 2058; 2059; 2060; 2061; 2062; 2063; 2064; 2065; 2066; 2067; 2068; 2069; 2070; 2071; 2072; 2073; 2074; 2075; 2076; 2077; 2078; 2079; 2080; 2081; 2082; 2083; 2084; 2085; 2086; 2087; 2088; 2089; 2090; 2091; 2092; 2093; 2094; 2095; 2096; 2097; 2098; 2099; 2100; 2101; 2102; 2103; 2104; 2105; 2106; 2107; 2108; 2109; 2110; 2111; 2112; 2113; 2114; 2115; 2116; 2117; 2118; 2119; 2120; 2121; 2122; 2123; 2124; 2125; 2126; 2127; 2128; 2129; 2130; 2131; 2132; 2133; 2134; 2135; 2136; 2137; 2138; 2139; 2140; 2141; 2142; 2143; 2144; 2145; 2146; 2147; 2148; 2149; 2150; 2151; 2152; 2153; 2154; 2155; 2156; 2157; 2158; 2159; 2160; 2161; 2162; 2163; 2164; 2165; 2166; 2167; 2168; 2169; 2170; 2171; 2172; 2173; 2174; 2175; 2176; 2177; 2178; 2179; 2180; 2181; 2182; 2183; 2184; 2185; 2186; 2187; 2188; 2189; 2190; 2191; 2192; 2193; 2194; 2195; 2196; 2197; 2198; 2199; 2200; 2201; 2202; 2203; 2204; 2205; 2206; 2207; 2208; 2209; 2210; 2211; 2212; 2213; 2214; 2215; 2216; 2217; 2218; 2219; 2220; 2221; 2222; 2223; 2224; 2225; 2226; 2227; 2228; 2229; 2230; 2231; 2232; 2233; 2234; 2235; 2236; 2237; 2238; 2239; 2240; 2241; 2242; 2243; 2244; 2245; 2246; 2247; 2248; 2249; 2250; 2251; 2252; 2253; 2254; 2255; 2256; 2257; 2258; 2259; 2260; 2261; 2262; 2263; 2264; 2265; 2266; 2267; 2268; 2269; 2270; 2271; 2272; 2273; 2274; 2275; 2276; 2277; 2278; 2279; 2280; 2281; 2282; 2283; 2284; 2285; 2286; 2287; 2288; 2289; 2290; 2291; 2292; 2293; 2294; 2295; 2296; 2297; 2298; 2299; 2300; 2301; 2302; 2303; 2304; 2305; 2306; 2307; 2308; 2309; 2310; 2311; 2312; 2313; 2314; 2315; 2316; 2317; 2318; 2319; 2320; 2321; 2322; 2323; 2324; 2325; 2326; 2327; 2328; 2329; 2330; 2331; 2332; 2333; 2334; 2335; 2336; 2337; 2338; 2339; 2340; 2341; 2342; 2343; 2344; 2345; 2346; 2347; 2348; 2349; 2350; 2351; 2352; 2353; 2354; 2355; 2356; 2357; 2358; 2359; 2360; 2361; 2362; 2363; 2364; 2365; 2366; 2367; 2368; 2369; 2370; 2371; 2372; 2373; 2374; 2375; 2376; 2377; 2378; 2379; 2380; 2381; 2382; 2383; 2384; 2385; 2386; 2387; 2388; 2389; 2390; 2391; 2392; 2393; 2394; 2395; 2396; 2397; 2398; 2399; 2400; 2401; 2402; 2403; 2404; 2405; 2406; 2407; 2408; 2409; 2410; 2411; 2412; 2413; 2414; 2415; 2416; 2417; 2418; 2419; 2420; 2421; 2422; 2423; 2424; 2425; 2426; 2427; 2428; 2429; 2430; 2431; 2432; 2433; 2434; 2435; 2436; 2437; 2438; 2439; 2440; 2441; 2442; 2443; 2444; 2445; 2446; 2447; 2448; 2449; 2450; 2451; 2452; 2453; 2454; 2455; 2456; 2457; 2458; 2459; 2460; 2461; 2462; 2463; 2464; 2465; 2466; 2467; 2468; 2469; 2470; 2471; 2472; 2473; 2474; 2475; 2476; 2477; 2478; 2479; 2480; 2481; 2482; 2483; 2484; 2485; 2486; 2487; 2488; 2489; 2490; 2491; 2492; 2493; 2494; 2495; 2496; 2497; 2498; 2499; 2500; 2501; 2502; 2503; 2504; 2505; 2506; 2507; 2508; 2509; 2510; 2511; 2512; 2513; 2514; 2515; 2516; 2517; 2518; 2519; 2520; 2521;
#aoristic weights for all sites with rainfall below 300mm
aoristicweights_below300mm <- sitecounts_below300mm <- below300mm
sitecounts_below300mm[newcolumns] <- NA
sitecounts_below300mm<-sitecounts_below300mm[,newcolumns]
aoristicweights_below300mm[newcolumns]<-NA
aoristicweights_below300mm<-aoristicweights_below300mm[,newcolumns]
aoristicweights_below300mm<- sitecounts_below300mm
for (i in 1:nrow(below300mm)){
cat(paste(i,"; ",sep=""))
sitestart <- below300mm$StartdatesBP[i]
siteend <- below300mm$EnddatesBP[i]
siteyears <- seq(siteend,sitestart,by=1)
siteyearshist <- hist(siteyears, breaks=breaks, plot=FALSE)
timeweights <- siteyearshist$counts/binwidth
sitecounts_below300mm[i,] <- timeweights #counts
aoristicweights_below300mm[i,] <- timeweights/sum(timeweights) #aoristic weights
}
## 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80; 81; 82; 83; 84; 85; 86; 87; 88; 89; 90; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100; 101; 102; 103; 104; 105; 106; 107; 108; 109; 110; 111; 112; 113; 114; 115; 116; 117; 118; 119; 120; 121; 122; 123; 124; 125; 126; 127; 128; 129; 130; 131; 132; 133; 134; 135; 136; 137; 138; 139; 140; 141; 142; 143; 144; 145; 146; 147; 148; 149; 150; 151; 152; 153; 154; 155; 156; 157; 158; 159; 160; 161; 162; 163; 164; 165; 166; 167; 168; 169; 170; 171; 172; 173; 174; 175; 176; 177; 178; 179; 180; 181; 182; 183; 184; 185; 186; 187; 188; 189; 190; 191; 192; 193; 194; 195; 196; 197; 198; 199; 200; 201; 202; 203; 204; 205; 206; 207; 208; 209; 210; 211; 212; 213; 214; 215; 216; 217; 218; 219; 220; 221; 222; 223; 224; 225; 226; 227; 228; 229; 230; 231; 232; 233; 234; 235; 236; 237; 238; 239; 240; 241; 242; 243; 244; 245; 246; 247; 248; 249; 250; 251; 252; 253; 254; 255; 256; 257; 258; 259; 260; 261; 262;
Combine and summarize the data into a format that can be graphed, colSum() will change the data.frame to numeric variables - this is done for both rainfall groups and all data.
#summaries of groups
sum_sitecounts_above300mm <- colSums(sitecounts_above300mm)
sum_aoristicweights_above300mm <- colSums(aoristicweights_above300mm)
sum_sitecounts_below300mm <- colSums(sitecounts_below300mm)
sum_aoristicweights_below300mm <- colSums(aoristicweights_below300mm)
sum_sitecounts<- colSums(sitecounts)
sum_aoristicweights <- colSums(aoristicweights)
Then need to extract the dates of use and convert into a format that can be compared.
#above 300mm (counts and weights)
sum_above300mm <- sum_sitecounts_above300mm[12:41] #extract values between 3000-6000BP, count
x <- sum_above300mm
sum_above300mm_norm<-((x-min(x))/(max(x) - min(x)))
sum_above300mm_aoristicwts <- sum_aoristicweights_above300mm[12:41] #weights
x <- sum_above300mm_aoristicwts
sum_above300mm_aoristicwts_norm<-((x-min(x))/(max(x) - min(x)))
#below 300mm(counts and weights)
sum_below300mm <- sum_sitecounts_below300mm[12:41] #counts
x <- sum_below300mm
sum_below300mm_norm<-((x-min(x))/(max(x) - min(x)))
sum_below300mm_aoristicwts <- sum_aoristicweights_below300mm[12:41] #weights
x <- sum_below300mm_aoristicwts
sum_below300mm_aoristicwts_norm<-((x-min(x))/(max(x) - min(x)))
#figure 4a (Raw count)
plot(mids[12:41],sum_above300mm_norm, lty = "solid", col = "white", cex.axis = 0.5, xlab = "cal BP", ylab = "", xlim = c(5900,3100))
lines(mids[12:41],sum_above300mm_norm, lty="solid", col="green")
lines(mids[12:41],sum_below300mm_norm, lty="solid", col="red")
rect(4300, -1, 3900, 1.5, col = rgb(128,128,128, alpha = 120, maxColorValue = 255), border = NA)
text(x = 4100, 0.90, labels = "4.2 k", font = 2, col = "black")
abline(v=seq(5800, 3200, -200), lty = 3, col= "grey86")
text(x = 5860, y = 0.6, labels = "Raw count", font = 2, cex = 0.5)
text(x = 5880, y = 0.55, label = paste("site-phase=", nrow(above300mm)," ", "sites=", length(unique(above300mm$Id)), sep = "", " (settlement)"), font = 1, col = "green", cex = 0.5, adj = c(0.1))
text(x = 5880, y = 0.50, label = paste("site-phase=", nrow(below300mm)," ", "sites=", length(unique(below300mm$Id)), sep = "", " (settlement)"), font = 1, col = "red", cex = 0.5, adj = c(0.1))
legend(x=5960,0.47, legend=c("Above 300mm","Sbelow 300mm"), lty=c("solid","solid"), lwd=c(0.5,0.5), col=c("green", "red"), bty="n", cex=0.5) #key
#figure 4b (Aoristic wts)
plot(mids[12:41], sum_above300mm_aoristicwts_norm, lty = "solid", col = "white", cex.axis = 0.5, xlab = "cal BP", ylab = "", xlim = c(5900,3100))
lines(mids[12:41], sum_above300mm_aoristicwts_norm, lty = "dashed", col = "green")
lines(mids[12:41], sum_below300mm_aoristicwts_norm, lty = "dashed", col = "red")
rect(4300, -1, 3900, 1.5, col = rgb(128,128,128, alpha = 120, maxColorValue = 255), border = NA)
text(x = 4100, 0.90, labels = "4.2 k", font = 2, col = "black")
abline(v=seq(5800, 3200, -200), lty = 3, col= "grey86")
text(x = 5800, y = 0.6, labels = "Aoristic Weights", font = 2, cex = 0.5)
text(x = 5880, y = 0.55, label = paste("site-phase=", nrow(above300mm)," ", "sites=", length(unique(above300mm$Id)), sep = "", " (settlement)"), font = 1, col = "green", cex = 0.5, adj = c(0.1))
text(x = 5880, y = 0.50, label = paste("site-phase=", nrow(below300mm)," ", "sites=", length(unique(below300mm$Id)), sep = "", " (settlement)"), font = 1, col = "red", cex = 0.5, adj = c(0.1))
legend(x=5960,0.47, legend=c("Above 300mm","Below 300mm"), lty=c("solid","solid"), lwd=c(0.5,0.5), col=c("green", "red"), bty="n", cex=0.5) #key
Figure 4 from paper (for comparison)
This dataset will be processed the same way as the other settlement dataset.
urban_sites <- read.csv("urban_sites.csv", sep = ",", header = T)
head(urban_sites)
## Name Id Rainfall StartDate EndDate SizeHa Longitude
## 1 Abu Harba 195 186.243 2500 1600 10.60000 37.38201
## 2 Abu Kula 190 437.290 2550 1500 10.00000 42.32779
## 3 Abu Kula 190 437.290 3100 2550 8.00000 42.32779
## 4 Agilah 221 255.632 1900 1650 31.87901 39.63867
## 5 Ain Dara 68 453.873 1600 1000 22.00000 36.85222
## 6 Al-Hirbat al-Garbiyya 193 185.758 1800 1600 22.00000 37.35318
## Latitude
## 1 35.44246
## 2 36.66782
## 3 36.66782
## 4 36.59434
## 5 36.46018
## 6 35.25884
## Source
## 1 PalaeoSyr Database (unbpublished but accessible to authors, see http://www.cepam.cnrs.fr/spip.php?rubrique175&lang=fr)
## 2 Wilkinson. T. J., Tucker. D. J., (1995) Settlement and Development in the North Jazira, Iraq: A study of the Archaeological Landscape, British Schhol of Archaeology in Iraq Department of Antiquities and Heritage, Baghdad.
## 3 Wilkinson. T. J., Tucker. D. J., (1995) Settlement and Development in the North Jazira, Iraq: A study of the Archaeological Landscape, British Schhol of Archaeology in Iraq Department of Antiquities and Heritage, Baghdad.
## 4 Corona Spy Photography Mission 1038, acquired 22/01/1967
## 5 Stone. E.C., Zimansky. P. E. (1999) The Iron Age Settlement at 'Ain Dara, Syria: Survey and Soundings, BAR International Series 786
## 6 PalaeoSyr Database (unbpublished but accessible to authors, see http://www.cepam.cnrs.fr/spip.php?rubrique175&lang=fr)
nrow(urban_sites) #283 occupation phases (corresponds with paper)
## [1] 283
#converted from BC/AD to BP
StartdatesBP_urban <- (urban_sites$StartDate) + 1950
EnddatesBP_urban <- (urban_sites$EndDate) +1950
urban_sites[, "StartdatesBP"] <- StartdatesBP_urban
urban_sites[, "EnddatesBP"] <- EnddatesBP_urban
#duration
Duration_urban <- urban_sites$StartdatesBP - urban_sites$EnddatesBP
urban_sites[, "Duration"] <- Duration_urban
head(urban_sites)
## Name Id Rainfall StartDate EndDate SizeHa Longitude
## 1 Abu Harba 195 186.243 2500 1600 10.60000 37.38201
## 2 Abu Kula 190 437.290 2550 1500 10.00000 42.32779
## 3 Abu Kula 190 437.290 3100 2550 8.00000 42.32779
## 4 Agilah 221 255.632 1900 1650 31.87901 39.63867
## 5 Ain Dara 68 453.873 1600 1000 22.00000 36.85222
## 6 Al-Hirbat al-Garbiyya 193 185.758 1800 1600 22.00000 37.35318
## Latitude
## 1 35.44246
## 2 36.66782
## 3 36.66782
## 4 36.59434
## 5 36.46018
## 6 35.25884
## Source
## 1 PalaeoSyr Database (unbpublished but accessible to authors, see http://www.cepam.cnrs.fr/spip.php?rubrique175&lang=fr)
## 2 Wilkinson. T. J., Tucker. D. J., (1995) Settlement and Development in the North Jazira, Iraq: A study of the Archaeological Landscape, British Schhol of Archaeology in Iraq Department of Antiquities and Heritage, Baghdad.
## 3 Wilkinson. T. J., Tucker. D. J., (1995) Settlement and Development in the North Jazira, Iraq: A study of the Archaeological Landscape, British Schhol of Archaeology in Iraq Department of Antiquities and Heritage, Baghdad.
## 4 Corona Spy Photography Mission 1038, acquired 22/01/1967
## 5 Stone. E.C., Zimansky. P. E. (1999) The Iron Age Settlement at 'Ain Dara, Syria: Survey and Soundings, BAR International Series 786
## 6 PalaeoSyr Database (unbpublished but accessible to authors, see http://www.cepam.cnrs.fr/spip.php?rubrique175&lang=fr)
## StartdatesBP EnddatesBP Duration
## 1 4450 3550 900
## 2 4500 3450 1050
## 3 5050 4500 550
## 4 3850 3600 250
## 5 3550 2950 600
## 6 3750 3550 200
#split into rainfall groups
urban_above300mm <- subset(urban_sites, urban_sites$Rainfall > 300)
urban_below300mm <- subset(urban_sites, urban_sites$Rainfall <= 300)
#basic descriptive stats of the rainfall groups
summary(urban_above300mm$Rainfall) #rainfall stats
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 304.2 350.4 422.9 425.5 474.4 735.9
summary(urban_below300mm$Rainfall)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 120.7 228.4 253.2 242.6 267.4 299.1
summary(urban_above300mm$Duration) #duration of site occupation in both groups
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 50.0 200.0 500.0 660.6 900.0 3000.0
summary(urban_below300mm$Duration)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 50 350 550 730 1000 3000
#bin into 100-year time slots
#reused same variables as above (binwidth, breaks, mids)
newcolumns_urban <- paste(breaks[1:(length(breaks)-1)], "start", sep = "")
#aoristic weights for urban sites (all)
#create df for results
aoristicweights_urban <- sitecounts_urban <- urban_sites
sitecounts_urban[newcolumns_urban] <- NA
sitecounts_urban<-sitecounts_urban[,newcolumns_urban]
aoristicweights_urban[newcolumns_urban]<-NA
aoristicweights_urban<-aoristicweights_urban[,newcolumns_urban]
aoristicweights_urban<- sitecounts_urban
for (i in 1:nrow(urban_sites)){
cat(paste(i,"; ",sep=""))
sitestart <- urban_sites$StartdatesBP[i]
siteend <- urban_sites$EnddatesBP[i]
siteyears <- seq(siteend, sitestart, by = 1)
siteyearshist <- hist(siteyears, breaks=breaks, plot = F)
timeweights <- siteyearshist$counts/binwidth
sitecounts_urban[i,] <- timeweights #counts
aoristicweights_urban[i,] <- timeweights/sum(timeweights) #aoristic weights
}
## 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80; 81; 82; 83; 84; 85; 86; 87; 88; 89; 90; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100; 101; 102; 103; 104; 105; 106; 107; 108; 109; 110; 111; 112; 113; 114; 115; 116; 117; 118; 119; 120; 121; 122; 123; 124; 125; 126; 127; 128; 129; 130; 131; 132; 133; 134; 135; 136; 137; 138; 139; 140; 141; 142; 143; 144; 145; 146; 147; 148; 149; 150; 151; 152; 153; 154; 155; 156; 157; 158; 159; 160; 161; 162; 163; 164; 165; 166; 167; 168; 169; 170; 171; 172; 173; 174; 175; 176; 177; 178; 179; 180; 181; 182; 183; 184; 185; 186; 187; 188; 189; 190; 191; 192; 193; 194; 195; 196; 197; 198; 199; 200; 201; 202; 203; 204; 205; 206; 207; 208; 209; 210; 211; 212; 213; 214; 215; 216; 217; 218; 219; 220; 221; 222; 223; 224; 225; 226; 227; 228; 229; 230; 231; 232; 233; 234; 235; 236; 237; 238; 239; 240; 241; 242; 243; 244; 245; 246; 247; 248; 249; 250; 251; 252; 253; 254; 255; 256; 257; 258; 259; 260; 261; 262; 263; 264; 265; 266; 267; 268; 269; 270; 271; 272; 273; 274; 275; 276; 277; 278; 279; 280; 281; 282; 283;
#aoristic weights for all sites with rainfall above 300mm
aoristicweights_urban_above300mm <- sitecounts_urban_above300mm <- urban_above300mm
sitecounts_urban_above300mm[newcolumns_urban ] <- NA
sitecounts_urban_above300mm <-sitecounts_urban_above300mm[,newcolumns_urban ]
aoristicweights_urban_above300mm[newcolumns_urban]<-NA
aoristicweights_urban_above300mm<-aoristicweights_urban_above300mm[,newcolumns_urban ]
aoristicweights_urban_above300mm<- sitecounts_urban_above300mm
for (i in 1:nrow(urban_above300mm)){
cat(paste(i,"; ",sep=""))
sitestart <- urban_above300mm$StartdatesBP[i]
siteend <- urban_above300mm$EnddatesBP[i]
siteyears <- seq(siteend,sitestart,by=1)
siteyearshist <- hist(siteyears, breaks=breaks, plot =F)
timeweights <- siteyearshist$counts/binwidth
sitecounts_urban_above300mm[i,] <- timeweights #counts
aoristicweights_urban_above300mm[i,] <- timeweights/sum(timeweights) #aoristic weights
}
## 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80; 81; 82; 83; 84; 85; 86; 87; 88; 89; 90; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100; 101; 102; 103; 104; 105; 106; 107; 108; 109; 110; 111; 112; 113; 114; 115; 116; 117; 118; 119; 120; 121; 122; 123; 124; 125; 126; 127; 128; 129; 130; 131; 132; 133; 134; 135; 136; 137; 138; 139; 140; 141; 142; 143; 144; 145; 146; 147; 148; 149; 150; 151; 152; 153; 154; 155; 156; 157; 158; 159; 160; 161; 162; 163; 164; 165; 166; 167; 168; 169; 170; 171; 172; 173; 174; 175; 176; 177; 178; 179; 180; 181; 182; 183; 184; 185; 186; 187; 188; 189; 190; 191; 192; 193; 194; 195; 196; 197; 198; 199; 200; 201; 202; 203;
#aoristic weights for all sites with rainfall below 300mm
aoristicweights_urban_below300mm <- sitecounts_urban_below300mm <- urban_below300mm
sitecounts_urban_below300mm[newcolumns_urban ] <- NA
sitecounts_urban_below300mm <-sitecounts_urban_below300mm[,newcolumns_urban ]
aoristicweights_urban_below300mm[newcolumns_urban]<-NA
aoristicweights_urban_below300mm<-aoristicweights_urban_below300mm[,newcolumns_urban ]
aoristicweights_urban_below300mm<- sitecounts_urban_below300mm
for (i in 1:nrow(urban_below300mm)){
cat(paste(i,"; ",sep=""))
sitestart <- urban_below300mm$StartdatesBP[i]
siteend <- urban_below300mm$EnddatesBP[i]
siteyears <- seq(siteend,sitestart,by=1)
siteyearshist <- hist(siteyears, breaks=breaks, plot =F)
timeweights <- siteyearshist$counts/binwidth
sitecounts_urban_below300mm[i,] <- timeweights #counts
aoristicweights_urban_below300mm[i,] <- timeweights/sum(timeweights) #aoristic weights
}
## 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58; 59; 60; 61; 62; 63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80;
#summaries of rainfall groups
sum_urbancounts <- colSums(sitecounts_urban)
sum_urbanaoristicweights <- colSums(aoristicweights_urban)
sum_urbancounts_above300mm <- colSums(sitecounts_urban_above300mm)
sum_urbanaoristicweights_above300mm <- colSums(aoristicweights_urban_above300mm)
sum_urbancounts_below300mm <- colSums(sitecounts_urban_below300mm)
sum_aoristicweights_below300mm <- colSums(aoristicweights_urban_below300mm)
#all urban sites (counts and weights converted)
urbancounts <- sum_urbancounts[12:41]
x <- urbancounts
urbancounts_norm <- ((x-min(x))/(max(x) - min(x)))
urbanwieghts <- sum_urbanaoristicweights[12:41]
x <- urbanwieghts
urbanwieghts_norm <- ((x-min(x))/(max(x) - min(x)))
#above 300mm (counts and weights converted)
urbansum_above300mm <- sum_urbancounts_above300mm[12:41] #extract values between 3000-6000BP
x <- urbansum_above300mm
urbansum_above300mm_norm<-((x-min(x))/(max(x) - min(x)))
urbansum_above300mm_aoristicwts <- sum_urbanaoristicweights_above300mm[12:41] #weights
x <- urbansum_above300mm_aoristicwts
urbansum_above300mm_aoristicwts_norm<-((x-min(x))/(max(x) - min(x)))
#below 300mm(counts and weights converted)
urbansum_below300mm <- sum_urbancounts_below300mm[12:41] #counts
x <- urbansum_below300mm
urbansum_below300mm_norm<-((x-min(x))/(max(x) - min(x)))
urbansum_below300mm_aoristicwts <-sum_aoristicweights_below300mm[12:41] #weights
x <- urbansum_below300mm_aoristicwts
urbansum_below300mm_aoristicwts_norm<-((x-min(x))/(max(x) - min(x)))
##Replication of Figure 5a
plot(mids[12:41],urbansum_above300mm_norm, lty = "solid", col = "white", cex.axis = 0.5, xlab = "cal BP", ylab = "", xlim = c(5900,3100))
lines(mids[12:41],urbansum_above300mm_norm, lty="solid", col="green")
lines(mids[12:41],urbansum_below300mm_norm, lty="solid", col="red")
rect(4300, -1, 3900, 1.5, col = rgb(128,128,128, alpha = 120, maxColorValue = 255), border = NA)
text(x = 4100, 0.4, labels = "4.2 k", font = 2, col = "black")
abline(v=seq(5800, 3200, -200), lty = 3, col= "grey86")
text(x = 5860, y = 0.6, labels = "Raw count", font = 2, cex = 0.5)
text(x = 5880, y = 0.55, label = paste("site-phase=", nrow(urban_above300mm)," ", "sites=", length(unique(urban_above300mm$Id)), sep = "", " (settlement)"), font = 1, col = "green", cex = 0.5, adj = c(0.1))
text(x = 5880, y = 0.50, label = paste("site-phase=", nrow(urban_below300mm)," ", "sites=", length(unique(urban_below300mm$Id)), sep = "", " (settlement)"), font = 1, col = "red", cex = 0.5, adj = c(0.1))
legend(x=5960,0.47, legend=c("Above 300mm","Below 300mm"), lty=c("solid","solid"), lwd=c(0.5,0.5), col=c("green", "red"), bty="n", cex=0.5) #key
#figure 5c (Aoristic wts)
plot(mids[12:41], urbansum_above300mm_aoristicwts_norm, lty = "solid", col = "white", cex.axis = 0.5, xlab = "cal BP", ylab = "", xlim = c(5900,3100))
lines(mids[12:41], urbansum_above300mm_aoristicwts_norm, lty = "dashed", col = "green")
lines(mids[12:41], urbansum_below300mm_aoristicwts_norm, lty = "dashed", col = "red")
rect(4300, -1, 3900, 1.5, col = rgb(128,128,128, alpha = 120, maxColorValue = 255), border = NA)
text(x = 4100, 0.90, labels = "4.2 k", font = 2, col = "black")
abline(v=seq(5800, 3200, -200), lty = 3, col= "grey86")
text(x = 5800, y = 0.6, labels = "Aoristic Weights", font = 2, cex = 0.5)
text(x = 5880, y = 0.55, label = paste("site-phase=", nrow(urban_above300mm)," ", "sites=", length(unique(urban_above300mm$Id)), sep = "", " (settlement)"), font = 1, col = "green", cex = 0.5, adj = c(0.1))
text(x = 5880, y = 0.50, label = paste("site-phase=", nrow(urban_below300mm)," ", "sites=", length(unique(urban_below300mm$Id)), sep = "", " (settlement)"), font = 1, col = "red", cex = 0.5, adj = c(0.1))
legend(x=5960,0.47, legend=c("Above 300mm","Below 300mm"), lty=c("solid","solid"), lwd=c(0.5,0.5), col=c("green", "red"), bty="n", cex=0.5) #key
Figure 5 from paper (for comparison)
First we load in the radiocarbon dataset and subset all values associated with short-lived dates as described in the paper.
library(rcarbon)
library(parallel)
#Load in radiocarbon data set and view first few rows
dates <- read.csv("dates.csv", sep = ",", header = T)
head(dates)
## DateID LabID OthLabID Problems CRA Error DC13 Material
## 1 6771 OxA-882 6100 120 0 grain (charred)
## 2 5075 Ly-12506 4070 40 0 charcoal
## 3 5076 Ly-12507 4090 35 0 seed
## 4 5077 Ly-12508 3990 40 0 seed
## 5 5078 Ly-12509 3990 40 0 seed
## 6 5079 Ly-12510 4020 40 0 seed
## Species SiteID SiteName
## 1 Triticum boeoticum 1060 Abu Hureyra
## 2 6606 Al Rawda
## 3 6606 Al Rawda
## 4 6606 Al Rawda
## 5 6606 Al Rawda
## 6 6606 Al Rawda
## SiteContext
## 1 Level 324 -basal pit-; Trench E, phase 1, basal pit; Trench E, level 468, pit
## 2 LOCUS 1084 - SECT. 1
## 3 LOCUS 1115 - SECT. 1
## 4 LOCUS 2036 - SECT. 2 A
## 5 LOCUS 2049 - SECT. 2 A
## 6 LOCUS 2063 - SECT. 2 A
## SiteType Country Longitude Latitude UTM_X UTM_Y LocQual
## 1 settlement SY 38.3933 35.8683 445228.5 3969511 B
## 2 settlement SY 37.6330 35.1810 375530.3 3893971 A
## 3 settlement SY 37.6330 35.1810 375530.3 3893971 A
## 4 settlement SY 37.6330 35.1810 375530.3 3893971 A
## 5 settlement SY 37.6330 35.1810 375530.3 3893971 A
## 6 settlement SY 37.6330 35.1810 375530.3 3893971 A
## Source Comments
## 1 Moore et al. 1986 ;CalPal; Flohr et al. 2016
## 2 BANADORA
## 3 BANADORA
## 4 BANADORA
## 5 BANADORA
## 6 BANADORA
#subset short lived dates (defined as a separate group in text)
shortlived_dates <- dates[grep(paste(c("antler","bone","coprolite","seed","fruit","carpal","tooth", "grain", "twig", "grape","monocot","collagen", "dung", "teeth","tissue", "plant remains","stone","onion"),collapse="|"), dates$Material), ]
Create parameters to work with for the analyses, including creating a available for bins of 50 years and creating a broader time span to avoid edge effects (as discussed in the paper).
detectCores()#detect number of cores in usage
## [1] 8
#Parameters for analyses
bin <- 50
runm <- 50 #for smoothing (number stipulated in text)
ncores <- 8
realstartBP <- 6000 #original time span as stipulated in paper
realendBP <- 3000
workingstartBP <- realstartBP + 500 #broader time span to avoid edge effect
workingendBP <- realendBP - 500
Then need to calibrate the radiocarbon dates using the calibrate() function from the ‘rcarbon’ package. The dates were calibrated as unnormalised as discussed previously. Although normalization (or not) does not have an impact on the shape of each individual dates calibrated probability distribution, it does influence the shape of SPDs. Additionally, a ‘IntCal20’ calibration curves was used as direct in the paper.
#Calibrate dates (un-normalised)
cal_dates <- calibrate(x = dates$CRA, errors = dates$Error, calCurves = "intcal20", normalised = F, ncores = ncores, calMatrix = F)
## [1] "Calibrating radiocarbon ages..."
## [1] "Running in parallel (standard calibration only) on 1 workers..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done."
cal_shortlived_dates <- calibrate(x = shortlived_dates$CRA, errors = shortlived_dates$Error, calCurves = "intcal20", normalised = F, ncores = ncores, calMatrix = F)
## [1] "Calibrating radiocarbon ages..."
## [1] "Running in parallel (standard calibration only) on 8 workers..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done."
#summary of calibrated dates for both groups
summary(cal_dates[1:10])
## DateID MedianBP OneSigma_BP_1 OneSigma_BP_2 OneSigma_BP_3 OneSigma_BP_4
## 1 1 6973 7157 to 7094 7082 to 6849 6813 to 6802 NA to NA
## 2 2 4563 4786 to 4766 4616 to 4589 4587 to 4517 4477 to 4446
## 3 3 4598 4793 to 4762 4689 to 4682 4622 to 4523 4457 to 4457
## 4 4 4469 4518 to 4468 4447 to 4417 NA to NA NA to NA
## 5 5 4469 4518 to 4468 4447 to 4417 NA to NA NA to NA
## 6 6 4484 4522 to 4458 4456 to 4436 4432 to 4424 NA to NA
## 7 7 4165 4235 to 4195 4187 to 4142 4129 to 4092 NA to NA
## 8 8 4285 4400 to 4370 4355 to 4325 4300 to 4234 4197 to 4185
## 9 9 4285 4400 to 4370 4355 to 4325 4300 to 4234 4197 to 4185
## 10 10 4369 4422 to 4348 4333 to 4295 NA to NA NA to NA
## OneSigma_BP_5 TwoSigma_BP_1 TwoSigma_BP_2 TwoSigma_BP_3 TwoSigma_BP_4
## 1 NA to NA 7256 to 6724 6708 to 6676 NA to NA NA to NA
## 2 NA to NA 4803 to 4757 4698 to 4673 4649 to 4423 NA to NA
## 3 NA to NA 4813 to 4753 4713 to 4666 4654 to 4515 4481 to 4444
## 4 NA to NA 4571 to 4542 4534 to 4388 4372 to 4353 4328 to 4299
## 5 NA to NA 4571 to 4542 4534 to 4388 4372 to 4353 4328 to 4299
## 6 NA to NA 4782 to 4768 4614 to 4595 4583 to 4411 NA to NA
## 7 NA to NA 4347 to 4335 4293 to 4076 4041 to 3990 NA to NA
## 8 4165 to 4161 4409 to 4223 4206 to 4155 NA to NA NA to NA
## 9 4165 to 4161 4409 to 4223 4206 to 4155 NA to NA NA to NA
## 10 NA to NA 4514 to 4481 4444 to 4284 4276 to 4247 NA to NA
summary(cal_shortlived_dates[1:10])
## DateID MedianBP OneSigma_BP_1 OneSigma_BP_2 OneSigma_BP_3 OneSigma_BP_4
## 1 1 6973 7157 to 7094 7082 to 6849 6813 to 6802 NA to NA
## 2 2 4598 4793 to 4762 4689 to 4682 4622 to 4523 4457 to 4457
## 3 3 4469 4518 to 4468 4447 to 4417 NA to NA NA to NA
## 4 4 4469 4518 to 4468 4447 to 4417 NA to NA NA to NA
## 5 5 4484 4522 to 4458 4456 to 4436 4432 to 4424 NA to NA
## 6 6 4165 4235 to 4195 4187 to 4142 4129 to 4092 NA to NA
## 7 7 5454 5580 to 5503 5489 to 5444 5406 to 5326 NA to NA
## 8 8 6445 6537 to 6516 6502 to 6310 NA to NA NA to NA
## 9 9 7352 7426 to 7272 NA to NA NA to NA NA to NA
## 10 10 7378 7504 to 7267 NA to NA NA to NA NA to NA
## TwoSigma_BP_1 TwoSigma_BP_2 TwoSigma_BP_3 TwoSigma_BP_4
## 1 7256 to 6724 6708 to 6676 NA to NA NA to NA
## 2 4813 to 4753 4713 to 4666 4654 to 4515 4481 to 4444
## 3 4571 to 4542 4534 to 4388 4372 to 4353 4328 to 4299
## 4 4571 to 4542 4534 to 4388 4372 to 4353 4328 to 4299
## 5 4782 to 4768 4614 to 4595 4583 to 4411 NA to NA
## 6 4347 to 4335 4293 to 4076 4041 to 3990 NA to NA
## 7 5658 to 5277 5172 to 5132 5104 to 5060 NA to NA
## 8 6667 to 6280 NA to NA NA to NA NA to NA
## 9 7557 to 7544 7508 to 7237 7221 to 7166 NA to NA
## 10 7612 to 7157 7106 to 7075 7035 to 7035 NA to NA
Summed probability distributions were calculated next. The function spd() aggregates calibrated radiocarbon dates within a defined chronological range (2,500 - 6,500 BP). To mitigate biases that can results from strong inter-site variability in sample size, artificial bins are implemented using the binPrep() function.
Normalised and non-normalised SPDs for both groups of dates were created in order to compare in Fig. 2c.
#Summed Probability Distributions
#SPD as normalised and non-normalised for comparison
bins <- binPrep(sites = dates$SiteID, age = dates$CRA, h = bin) #create bins
spdnorm <- spd(x = cal_dates, bins = bins, timeRange = c(workingstartBP, workingendBP), datenormalised = T, runm = runm)
## [1] "Extracting and aggregating..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done."
spd <- spd(x = cal_dates, bins = bins, timeRange = c(workingstartBP, workingendBP), datenormalised = F, runm = runm)
## [1] "Extracting and aggregating..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done."
#short-lived dates SPD
bins_shortlived <- binPrep(sites = shortlived_dates$SiteID, age = shortlived_dates$CRA, h = bin)
spdnorm_shortlived <- spd(x = cal_shortlived_dates, bins = bins_shortlived, timeRange = c(workingstartBP, workingendBP), datenormalised = T, runm = runm)
## [1] "Extracting and aggregating..."
##
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## [1] "Done."
spd_shortlived <- spd(x = cal_shortlived_dates, bins = bins_shortlived, timeRange = c(workingstartBP, workingendBP), datenormalised = F, runm = runm)
## [1] "Extracting and aggregating..."
##
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## [1] "Done."
Pearsons coefficient is then calculated. The Pearsons coefficient matches that described in the paper, r = 0.96 for both date groups.
#Calculate pearsons coefficient, both correlated to 0.96 as per the paper
#Pearsons coefficient for SPD of all dates (normalised and not)
(cor.test(spd$grid$PrDens, spdnorm$grid$PrDens, method = "pearson"))
##
## Pearson's product-moment correlation
##
## data: spd$grid$PrDens and spdnorm$grid$PrDens
## t = 247.17, df = 3999, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9668321 0.9706431
## sample estimates:
## cor
## 0.9687948
#Pearsons coefficient for SPD of short-lived dates (normalised and not)
(cor.test(spd_shortlived$grid$PrDens, spdnorm_shortlived$grid$PrDens, method = "pearson"))
##
## Pearson's product-moment correlation
##
## data: spd_shortlived$grid$PrDens and spdnorm_shortlived$grid$PrDens
## t = 235.17, df = 3999, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9635404 0.9677234
## sample estimates:
## cor
## 0.9656945
#Visualizing Bins (figure 2b)
dates_bins_med <- binMed(x = cal_dates, bins = bins)
## [1] "Extracting and aggregating..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done"
hist(dates_bins_med, 50, main="Histogram of number of bins", xlab = "cal BP", xlim = c(6000, 3000))
#plot unnormalised and normalised SPD (Figure 2c)
plot(spd, xlim=c(6000,3000),rescale = T) #spd data
plot(spdnorm, xlim=c(6000,3000), xaxt='n', yaxt='n', type='simple', col="black", rescale = T, add = T)#spd normalised data
plot(spd_shortlived, xlim=c(6000,3000), xaxt='n', yaxt='n', type='simple', col="green", rescale = T, add = T) #spd short-lived (not normalised)
abline(v=seq(5800, 3200, -200), lty = 3, col= "white") #mark every 200 years with dotted line
text(x = 5900, 1, labels="SPDs of calibrated radiocarbon dates", font=2, cex=0.75, adj=c(0,1)) #subtitle
text(x = 5800, 0.90, labels = paste("n =", nrow(dates), "sites = ", length(unique(dates$SiteID)), "bins = ", length(unique(bins)), "(all dates)"), font = 1, cex = 0.70, adj = c(0.1)) #data info (all dates)
text(x = 5780, 0.85, labels = paste("n =", nrow(shortlived_dates), "sites = ", length(unique(shortlived_dates$SiteID)), "bins = ", length(unique(bins_shortlived)), "(shortlived dates)"), font = 1, cex = 0.70, adj = c(0.1)) #data info (short-lived dates)
legend(x=5950,0.80, legend=c("SPD unnormalised","SPD normalised","SPD unnormalised (shortlived)"), lty=c("solid","solid","solid"), lwd=c(3,0.5,0.5), col=c("grey","black","green"), bty="n", cex=0.8) #key
A logitsic growth model is then fitted to the observed SPD and a 95% confidence envelope of 1,000 SPDs is randomly generated. This will statistically test if the observed pattern differs from the null models. As discussed in the paper, deviations outside of the 95% confidence limits indicate periods of population growth or decline that are greater than expected according to a logistic model of population growth.
The text did not stipulate which function to use to get a logistic model. I ran the the genetic logistic model test using the function modelTest() but the results did not compare to the paper. I then tried a custom logistic model and was successful.
The p-value is the same as that detailed in the paper (<0.01).
#testing against theoretical growth model
#Generate Logistic model
set.seed(1)
nsim <- 1000 #from text
#logistic_model <- modelTest(cal_dates, errors = dates$Error, bins = bins, nsim = nsim, timeRange = c(6000,3000), model = "exponential", runm = runm)
#tried this but didn't match the papers results so did a custom Log mod
#Custom Logistic model
spd.smoothed <- spd(cal_dates, timeRange = c(6000, 3000), bins = bins, runm = runm) #generate a smoothed SPD
## [1] "Extracting and aggregating..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done."
logFit <- nls(PrDens~SSlogis(calBP, Asym, xmid, scale), data = spd.smoothed$grid, control = nls.control(maxiter = 200), start = list(Asym = 0.2, xmid = 5500, scale = -100))
#generate a df containing the fitted values
logFitDens=data.frame(calBP=spd.smoothed$grid$calBP,PrDens=SSlogis(input=spd.smoothed$grid$calBP,Asym=coefficients(logFit)[1],xmid=coefficients(logFit)[2],scal=coefficients(logFit)[3]))
# Use the modelTest function
Log_model <- modelTest(cal_dates, errors = dates$Error, bins = bins, nsim = nsim, timeRange = c(6000,3000), model = "custom", predgrid = logFitDens, runm = runm, ncores = ncores, raw = T)
## [1] "Aggregating observed dates..."
## [1] "Monte-Carlo test..."
## [1] "Running in parallel on 8 workers..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Done."
# Retrieve p-values
p <- Log_model$pval
#plot log model (figure 2d)
plot(Log_model, xlim=c(6000,3000))
lines(Log_model$fit$calBP, Log_model$fit$PrDens, col = "black", lty = 2) #regression line
abline(v=seq(5800, 3200, -200), lty = 3, col= "white") #mark every 200 years with dotted line
abline(v=c(4300, 3900), col = c("red", "red"), lwd = c(2, 2)) #4.2k marker
text(x = 4100, 0.26, labels = "4.2 k", font = 2, col = "red")
text(x = 5900, 0.27, labels="Logistic Fit", font=2, cex=0.75, adj=c(0,1)) #subtitle text
legend(x=5950,0.26, legend=c("SPD (dates not normalised)", "Logistic Model", "95% MC envelope", "positive deviation", "negative deviation"), lty=c("solid","dashed","solid", "solid", "solid"), lwd=c(0.5,0.5, 3, 3, 3), col=c("black", "black", "grey", "red", "blue"), cex=0.7, bty="n") #key
text(x = 5700, 0.16, cex = 0.6, labels = substitute(paste("p value = ", x), list(x = round(Log_model$pval,4)))) #insert p value
Figure 2 from paper (for comparison)
The replicated values are almost identical to those produced in the paper. There is a slight difference in the pearsons value between SPD when compared to the count and aoristic, although this number is only .01 of a position off in both cases and can still be interpreted in the same light as the paper.
x <- c(spd$grid[3501:502,2]) #group values (late to recent)
SPD <- colSums(matrix(x, nrow = 100, ncol = 30)) #create matrix
count_1 <- sum_sitecounts[12:41] #extract dates of value
aoristic_1 <- sum_aoristicweights[12:41] #extract dates of value
proxies_1 <- data.frame(cbind(count_1, aoristic_1, SPD)) #combine variables into df
(table_3_replication <- cor(proxies_1, method = "pearson")) #peasons correlation between variables
## count_1 aoristic_1 SPD
## count_1 1.00000000 0.03816917 -0.1705716
## aoristic_1 0.03816917 1.00000000 0.2636534
## SPD -0.17057162 0.26365336 1.0000000
Table 3 from paper (for comparison)
Again, the replicated values are almost identical to those produced in the paper. There is a slight difference in the pearsons value between SPD when compared to the count and aoritsic, although this difference is slight and the signfiaince of these values compare to those presented in the paper.
x <- c(spd$grid[2501:502,2])
SPD <- colSums(matrix(x, nrow = 100, ncol = 20))
count_2 <- sum_sitecounts[22:41]
aoristic_2 <- sum_aoristicweights[22:41]
proxies_2 <- data.frame(cbind(count_2, aoristic_2, SPD))
(table_4_replication <- cor(proxies_2, method = "pearson"))
## count_2 aoristic_2 SPD
## count_2 1.0000000 -0.3894488 -0.5911249
## aoristic_2 -0.3894488 1.0000000 0.9176673
## SPD -0.5911249 0.9176673 1.0000000
#check that p value for one of the bold numbers in table 4, i.e correlation between SPD and aoristic.
#the p-value is <0.05 and is significant
(cor.test(SPD, aoristic_2, method = "pearson"))
##
## Pearson's product-moment correlation
##
## data: SPD and aoristic_2
## t = 9.7983, df = 18, p-value = 1.221e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8000270 0.9673571
## sample estimates:
## cor
## 0.9176673
Table 4 from paper (for comparison)
The replicated values are almost identical to those produced in the paper. There is are very minute difference in the pearsons value between SPD when compared to the count and aoritsic, however this correlation can be interpreted in the same way as the paper.
x <- c(spd$grid[3501:502,2])
SPD <- colSums(matrix(x, nrow = 100, ncol = 30))
count_3 <- sum_urbancounts[12:41]
aoritsic_3 <- sum_urbanaoristicweights[12:41]
proxies_3 <- data.frame(cbind(count_3, aoritsic_3, SPD))
(Table_5_replication<- cor(proxies_3, method = "pearson"))
## count_3 aoritsic_3 SPD
## count_3 1.0000000 0.8770329 0.8555891
## aoritsic_3 0.8770329 1.0000000 0.6669123
## SPD 0.8555891 0.6669123 1.0000000
#check that p value for one of the bold numbers in table 5, i.e correlation between SPD and count.
#the p-value is <0.05 and is significant
(cor.test(SPD, count_3, method = "pearson"))
##
## Pearson's product-moment correlation
##
## data: SPD and count_3
## t = 8.7459, df = 28, p-value = 1.699e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7160309 0.9293830
## sample estimates:
## cor
## 0.8555891
Table 5 from paper (for comparison)
The replicated values are identical to the those presented in the paper.
x <- c(spd$grid[2501:502,2])
SPD <- colSums(matrix(x, nrow = 100, ncol = 20))
count_4 <- sum_urbancounts[22:41]
aoritsic_4 <- sum_urbanaoristicweights[22:41]
proxies_4 <- data.frame(cbind(count_4, aoritsic_4, SPD))
(Table_5_replication<- cor(proxies_4, method = "pearson"))
## count_4 aoritsic_4 SPD
## count_4 1.0000000 0.9767146 0.9490895
## aoritsic_4 0.9767146 1.0000000 0.9416295
## SPD 0.9490895 0.9416295 1.0000000
#check that p value for one of the bold numbers in table 6, i.e correlation between SPD and count.
#the p-value is <0.05 and is significant
(cor.test(SPD, count_4, method = "pearson"))
##
## Pearson's product-moment correlation
##
## data: SPD and count_4
## t = 12.783, df = 18, p-value = 1.813e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8733820 0.9800129
## sample estimates:
## cor
## 0.9490895
Table 6 from paper (for comparison)
Overall, the analyses that was replicated in this report was successful. Both archaeological datasets were successfully loaded in and converted into a format that was then used for analyses. Basic descriptive stats for duration and rainfall variables in each dataset were calculated. The data was successfully binned into the time slices outlined in the paper and aoristic weights were successfully calculated using a loop. The process of calculating these weights using r code was not detailed in the paper which made it challenging to replicate. However, the authors did explain the theory behind getting aoristic weights of each site, enough so, that creating a loop to calculate these values was achievable. Replication of parts of the archaeological dataset figures (Fig. 4 and 5) was fairly straight forward and although each figure consisted of a lot of elements, it was not challenging to replicate.
While the paper has sufficient information on processing of the radiocarbon dataset, it was still a challenge to replicate. To assist in this data processing, I found the rcarbon vignette extremely useful (https://cran.r-project.org/web/packages/rcarbon/vignettes/rcarbon.html) and as a result, replication was successful.
The main challenges I faced was trying to process the spatial datasets for analyses. I was able to load the various shp files using the {rgdal} package, however further processing of the files was unclear, and I was unsuccessful in my attempt. A further challenge related to this issue was converting the datasets into an object that could be spatial mapped. While I think the next step would be to extract the coordinates of each site to then spatially map, there was no instruction of how this could be achieved in r. The issue of spatial mapping came up in the radiocarbon part of the paper, in which the dates are divided into rainfall groups. The authors appear to group these dates using the rainfall_zone.shp file, however, no further detail was provided on the specific steps to achieve this and I was unsuccessful in this part of the replication.
Apart for these challenges, all other aspects of replication was successful. The only difference I could spot between the original and replicated results were the very slight differences in the pearson coefficient values (Table 3 and 4). However, I think these differences were not significance enough to consider a failed result. I speculate they might have occurred due to outliers that the original methodology might have excluded. Overall, the results I achieved seemed to align with those presented in the paper. Given more time (and knowledge on the subject) it would also be interesting to replicate the spatial analyses.
Lawrence D, Palmisano A, de Gruchy MW (2021) Collapse and continuity: A multi-proxy reconstruction of settlement organization and population trajectories in the Northern Fertile Crescent during the 4.2kya Rapid Climate Change event. PLoS ONE 16(1): e0244871. https://doi.org/10.1371/journal.pone.0244871